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The  purpose  of  this  study  was  to  develop  a  computer 
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Abstract 


The  purpose  oi  this  study  was  to  develop  a  computer 
program  that  incorporates  efficient  techniques  tor  solving 
optimization  of  experimental  test  or  simulation  models  The 
program  is  interactive  and  user-friendly.  The  program  is 
written  in  Fortran  but  can  be  attach  to  any  simulation  model 
or  experiment.  The  program  is  limited  to  two  independent 
variables  and  one  dependent  variable.  The  a  1 gor 1 thm  o 1  t he 
main  program  is  steepest  a s c e n t  p a r t a n . 

The  study  compared  several  gradient  methods  and  round 

2 - k  factorial  the  most  efficient.  The  study  also  concluded 

that  Davies,  Swann,  and  Campey  *  DSC .'  -  F  owe  1  1  was  the  most 

useful  line  search.  The  study  uses  an  improvement  by  Faix 

t:  the  FAriTAN  method  to  eliminate  the  !  ma.  ,;ne  search. 

T  h e  p r ogr am  is  den lg n  ed  to  e 1  I l c l e n  1 1 y  solve  s  e  c o n i- order 

equations  and  iess  efficiently  higher  order  equations. 
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SEARCH:  AN  INTERACTIVE  COMPUTER  PROGRAM 
FOR  OPTIMIZING  TWO- VARI ABLE ,  UNCONSTRAINED 


EXPERIMENTS  OR  SIMULATION 


In  science,  industry,  military  operations,  business  and 
most  facets  of  life,  countless  effort  is  spent  trying  to 
improve  "production"  and  maximize  one's  output.  For  some 
examples,  in  chemistry,  what  amount  of  chemical  A  mixed  with 
chemical  B  produces  the  largest  amount  of  product  C.  In 
military  operations,  how  many  and  of  what  type  weapon  should 
attack  which  targets  to  optimize  the  damage  expectancy.  In 
business,  how  many  of  a  certain  item  type  should  be  ordered 
to  minimize  storage  while  maximizing  sales.  Or,  in  one’s 
personal  finance,  what  amount  should  be  paid  on  certain 
debts  and  what  amounts  should  be  invested  in  which  tvpe  of 
savings  to  maximize  one’s  wealth.  So,  whether  it  is  the 
right  amount  of  an  input  that  produces  the  best  output  or 
the  correct  number  of  inputs  to  provide  the  best  results, 
optimization  is  a  problem  dealt  with  almost  daily.  Hillier 
and  Lieberman  point  out,  "this  'search  for  optimality'  is  a 
very  important  theme  in  operations  research"  (3:5). 

Operations  research  covers  a  broad  spectrum  of 


optimization  programming.  At  one  end  lies  the  well 
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developed  linear  programming  with  its  maximization  of  a 
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linear  objective  function  and  restrictions  by  linear 
constraints.  Going  down  the  spectrum  in  limitations,  one 
finds  the  lesser  developed  non-linear  programs  such  as 
quadratic,  geometric,  and  fractional  programming.  These 
algorithms  still  optimize  an  objective  function,  but  the 
restrictions  of  linearity  are  relaxed.  Still  further  down 
this  restrictive  ladder,  one  comes  to  direct  search  methods. 
Here,  optimization  of  the  problem  is  still  the  goal,  but  the 
formula  of  the  objective  function  is  unknown.  Direct  search 
techniques  use  experiments  (trial  and  error)  to  gain 
information  about  the  optimum.  This  method  involves  either 
experimenting  with  the  real  system  itself  or  experimenting 
with  a  simulation  model  of  the  system.  This  search  category 
includes  techniques  such  as  exhaustive  search,  random 
search,  and  direct  search. 

Wilde  clarifies  this  category  of  problems  as  follows: 


The  search  problem  is  to  find,  after  only  a  few 
experiments,  a  set  of  operating  conditions  yielding  a 
value  of  the  criterion  y  which  is  close  to  the  best 
attainable.  From  another  point  of  view,  the  problem  is 
to  reach  a  specified  minimum  acceptable  level  of 
performance  in  as  few  trials  as  possible 
Geometrically  speaking,  we  would  like  to  climb  up  the 
r<  sponse  surface  as  quickly  as  possible,  even  though 
the  only  information  we  have  about  the  surface  comes 
from  the  past  experiments  we  have  run  [9:64]. 


New  and  better  search  techniques  are  evolving  as  the 
research  for  more  efficient  methods  is  continued.  The 
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thrust  of  this  study  and  criterion  of  effectiveness  of 
methods  is  efficiency  in  reaching  the  optimum.  More 
efficient  is  defined  as  requiring  fewer  experiments  to  reach 
the  vicinity  of  the  optimum.  Thus,  better,  more  efficient 
methods  arrive  at  the  optimum  using  fewer  experiments  or 
s  1  mu  1  a  1 1  o  n  s  . 

The  simplest,  and  probably  first,  search  technique  used 
was  the  exhaustive  search,  or  mere  enumeration.  This 
technique  involves  looking  at  all  possible  combinations  of 
input  variables  and  selecting  the  combination  giving  the 
highest  output.  This  accurate  method  probably  saw  a  short 
rebirth  with  the  advancements  of  computers.  The  exhaustive 
method  works  fine  for  small  problems;  however,  even  with 
computers,  a  more  efficient  method  is  needed  to  save  time 
and  cost. 

Since  with  exhaustive  search  it  might  often  times  be 
prohibitive,  random  search  might  be  used.  This  involves 
randomly  selecting  input  combinations  for  testing.  The 
problem  with  random  search  is  never  knowing  when  the  optimum 
has  been  reached  unless  all  points  are  tested.  If  only  a 
few  experiments  are  possible,  random  search  might  be 
considered  the  best  choice  for  a  large  problem. 

From  the  need  for  a  more  scientific  method,  the  direct 
search  methods  were  developed  Direct  search  is  a  planned, 
mathematical  search  that  leads  one  to  the  optimum.  Over  the 
last  thirty  years  there  have  been  numerous  direct  search 


plans  developed.  Search  plans  basically  fall  into  tv;o 
categories:  simultaneous  and  sequential.  Plans  specifying 

the  location  of  every  experiment  before  any  results  are 
known  are  called  simultaneous,  while  a  plan  permitting  the 
experimenter  to  base  future  experiments  on  past  outcomes  are 
called  sequential  (8:5)  .  Simultaneous  search  plans,  usually 
called  experi me ntal  designs,  have  been  developed  that 
s  y  s  t  e  ma  tically  test  points  in  a  specified  region  of 
interest.  Response  surface  methodology  (RSM)  takes  the 
experimental  design  and  calculates  an  estimated  equation  of 
the  real  system  from  which  an  expected  optimum  can  be 
derived  mathematically.  Numerous  sequential  search 
algorithms  have  also  been  developed.  These  algorithms 
generally  entail  the  use  of  gradients,  directional  line 
searches,  geometry,  and  sometimes  curve  fitting. 

Problem  _S_t  at  emend, 

Frequently,  engineers  are  given  problems  to  solve  in 
which  they  want  the  optimal  solution  (either  maximization  or 
minimization)  and  no  equations  or  formulas  exist  of  the 
objective  function.  Thus  experimenting  (or  simulating) 
provides  the  only  clues  for  the  location  of  the  optimum. 

There  is  a  need  for  a  computer  software  program  that 
guides  a  person  to  the  optimum  whether  maximum  or  minimum 
when  only  simulation  or  experimentation  is  available.  It 
should  combine  various  efficient  direct  search  techniques  in 


an  algorithm  to  expedite  the  search.  It  should  employ 
simple  techniques  that  a  practicing  engineer  should 
understand.  It’s  goal  should  be  to  minimize  the  number  of 
experimental  tests  (simulations)  required.  Basically,  it 
should  be  simple  for  the  user  to  implement  and  use. 

Research.  Objectives 

The  overall  objective  of  this  research  is  to  develop  an 
interactive,  user-friendly  computer  package  which  allows  one 
to  find  the  optimal  response  to  experimental  test  or 
simulation  models.  The  program  will  contain  the  most 
efficient  search  techniques.  The  program  will  quickly  solve 
for  the  optimum  for  quadratic  surfaces  and  many  higher  order 
equat ions. 

Subobjectives  of  this  research  effort  are  as  follows: 

(1)  The  efficiency  of  different  techniques  during 
different  phases  of  the  program  compared  In  order  to  select 
the  most  efficient  techniques  for  the  program.  The  measure 
of  effectiveness  for  efficiency  being  the  least  experimental 
trials  needed  for  the  required  accuracy. 

(2)  A  verification  of  the  co mp uter  program 
accomplished  showing  that  it  does  solve  optimization 
problems.  This  would  entail  taking  various  problems  and 
checking  to  see  if  the  program  can  find  the  optimal 
solution  . 


The  limitations  of  this  research  are: 


(1)  The  independent  (input)  variables  are  limited  to 
two  and  they  are  continuous  real  variables. 

(2)  There  is  a  single  dependent  (output)  variable  from 
the  model  considered  in  the  program. 

(3)  The  only  experimental  designs  used  are  a  2-K 
factorial  for  first-order  equation  fit  and  a  3-K  factorial 
for  second-order  equation  fit. 

(4)  Experimental  error  is  mainly  handled  by  repeating 
the  simulation  test  and  then  averaging  the  results.  The 
number  of  repetitions  are  at  the  discretion  of  the  user 

(5)  The  validation  of  the  user  friendliness  of  the 
program  is  accomplished  by  having  a  fellow  student  run  the 
program  unaided  to  solve  a  problem. 


S  umma  r.y 

This  chapter  briefly  discussed  the  general  background, 
problem  state  me  nt,  research  objectives,  and  scope  pertaining 
to  this  research.  The  next  chapter  will  discuss  the 
literature  and  methodology  pertinent  to  this  research 


effort. 


II.  Literature  And  Methodology 


In  the  past,  the  practical  method  known  for  handling 
optimization  problems  was  the  classical  differential 
calculus.  However,  the  classical  method  is  impossible  to 
use  when  the  objective  function  is  undefined.  Therefore,  an 
indirect  method  of  finding  the  optimum  using  trial  and  error 
must  be  used.  Wilde  gives  the  name  "optimum  seeking 
procedures'  to  the  strategies  guiding  search  for  the  optimum 
of  any  function  about  which  full  knowledge  is  not  available 
(  9  :  v  1  i  1  )  . 

Plan  QL  Attack 

The  only  way  to  gain  information  about  an  unknown 
function  is  by  direct  measurement,  in  other  words 
experimentation  (5:vii) .  In  this  optimum  seeking  method, 
each  experiment  has  two  purposes,  not  only  to  attain  a  good 
response  surface  value,  but  also  to  give  information  useful 
for  locating  future  experiments  where  desirable  values  of 
the  response  surface  are  likely  to  be  found.  Thus, 
throughout  the  search  one  must  continually  be  deciding  to 
climb  or  to  explore.  At  the  beginning,  when  nothing  at  all 
is  known  about  the  function,  one  must  explore  in  some  small 
region,  usually  chosen  as  a  best  guess,  so  that  one  might 
place  the  following  experiments  in  an  uphill  direction.  In 
the  middle  of  the  search,  after  having  explored  some  region. 


one  tries  to  climb  as  fast  as  possible,  exploring  only  when 
strictly  necessary  to  guide  the  successive  steps.  Toward 
the  end  of  the  search,  when  one  is  finally  near  the  top, 
extensive  exploration  may  be  needed  to  attain  any  increase 
in  elevation,  the  slope  of  the  response  surface  often  being 
slight  near  the  optimum  (9:64). 

An  analogy  to  this  plan  of  attack  might  be  like  a  blind 
man  climbing  to  the  top  (highest  point)  of  a  mountain.  At 
the  bottom  of  the  mountain,  he  probes  around  with  his  cane 
to  find  the  steepest  uphill  slope.  After  this  initial 
exploring,  he  proceeds  in  this  uphill  direction  until  he 
reaches  a  point  where  he  starts  to  go  downhill.  At  this 
point  he  probes  around  this  area  for  a  new  uphill  direction 
and  proceeds  uphill  again.  This  continues  until  he  reaches 
the  top  and  can  find  no  new  uphill  direction.  At  this  time, 
he  explores  extensively  around  the  top  to  find  the  upmost 
highest  point.  The  direct  search  method  is  similar  to  this 
blind  man’s  search  in  that  one  cannot  see  where  one  is 
going,  but  only  by  probing  with  experiments,  like  searching 
with  a  cane,  can  one  get  the  direction  to  the  optimum. 

The  three  phases  of  the  attack  plan  will  now  be 
discussed  separately.  The  beginning  exploratory  phase  to 
find  the  uphill  direction  will  be  called  the  gradient  phase. 
The  middle  climbing  phase  will  be  called  the  acceleration 
phase.  The  final  phase  exploring  the  top  will  be  called  the 
second-order  exploratory  phase 


Gradi  ant  E h.as  e 


In  this  phase  of  the  search,  one  is  exploring  for  the 
uphill  gradient  (slope)  of  the  response  surface  from  the 
initial  starting  point.  As  mentioned  earlier,  if  the 
function  is  known,  derivatives  could  be  taken  at  this  point 
to  find  the  gradient.  Since  it  is  unknown,  another  way  must 
be  determined  to  attain  the  gradient.  Note--This  is  one  of 
the  mam  differences  between  some  non-linear  programs  and 
these  direct  search  techniques. 

Wilde  proposes  one  method  of  obtaining  the  general 
slope  of  the  response  surface  in  the  neighborhood  of  the 
initial  point.  First,  find  the  gradient  in  the  xl -direction 
parallel  to  the  xl  axis.  To  do  this,  one  varies  the  xl 
value  slightly  while  holding  the  x2  coordinate  at  x20  :  the 

initial  x2  value)  allowing  just  enough  distance  between  xll 
and  xlO  (the  initial  xl  value)  to  make  the  outcome  yl 
distinguishable  from  yO  (the  initial  y  response  v a  1 u e i  The 

straight  line  through  the  points  yO  and  yl  lies  entirely  in 
the  plane  of  the  x20  value  and  is  approximately  t anger.*  ♦.  : 
the  response  surface  at  yO .  The  slope  of  this  line  is  given 
by  yl-yO/xll-xlO.  Now,  a  similar  exploratory  experiment  is 
done,  but  this  time  varying  the  x2  coordinate  and  holding 
the  xl  coordinate  constant  at  xlO.  The  straight  line 
passing  through  y  2  and  yO  lies  entirely  in  the  plane  of  xlO 
and  the  slope  of  this  line  is  given  by  y2-v0/x22-x20.  The 


I  I  -  9 


three  points  yO ,  yl,  and  y2  on  the  response  surface  are 
enough  to  determine  the  plane  approximately  tangent  to  the 
surface  at  yO .  An  equation  of  this  plane  would  be  of  the 
form  y=  bO + b 1 x 1 + b2x2  (9:65-68) .  From  the  above,  one  can 
deduce  the  direction  to  proceed  from  the  initial  point  is  on 
a  vector  that  goes  y2-y0  in  the  x2  direction  while  going  yl- 
yO  in  the  xl  direct,. on.  Thus,  with  only  two  experiments  an 
approximate  direction  to  start  the  climb  is  found. 

Another  exploratory  search  method  which  uses  from  two 
to  four  experiments  is  given  by  R.  Hooke  and  T.A.  Jeeves. 
This  method  is  similar  to  the  previous  Wilde  method,  but  is 
more  sequential.  The  gradient  is  obtained  as  follows:  After 
the  initial  point  (xl0,x20)  is  evaluated  for  yO .  xlO  is 
changed  by  an  incremental  amount,  +  rf,  so  that  x 1  1 = x  1  0  * r f 
If  the  y  response  is  an  improvement  over  yO.  then  xll  is 
adopted  as  the  new  coordinate  in  the  xl  direction  If 
xlO+rf  fails  to  improve  the  response,  xlO  is  changed  bv  -rf 
and  the  value  of  the  y  response  again  checked  for 
improvement.  If  the  value  of  y  is  not  improved  by  either 
xlO  *  rf,  xlO  is  left  unchanged  After  the  xl  d  :  r  e  c  *  i  '  r.  :  s 
modified,  then  x20  is  changed  by  an  amount,  ♦  r  f  ,  and  the 
above  test  is  repeated  in  the  x2  direction  to  complete  :  r.  e 
exploratory  search.  The  successfully  changed  variables 
define  a  vector  from  the  initial  point  for  a  direction  to  do 


an  acceleration  phase  (  4  ■  1  4  2  -  1  4  8  ) 


1  0 


Another  method  to  find  the  gradient  is  to  use  a  2-k 
factorial  design  and  fit  a  first-order  equation  to  the 
response  surface  using  RSM  procedures.  What  are  the 
advantages  of  using  a  factorial  design7  Montgomery 
concludes  that  factorial  designs  are  more  efficient  than 
one - f ac tor -at  -  a- 1 1  me  experiments.  Also,  a  factorial  design 
is  necessary  when  interactions  may  be  present,  to  avoid 
misleading  conclusions.  Finally,  factorial  designs  allow 
effects  of  a  factor  to  be  estimated  at  several  levels  of  the 
other  factor,  yielding  conclusions  that  are  valid  over  a 
range  of  experimental  conditions  (5: 192) .  A  2-k  design 
would  be  like  a  box  drawn  around  the  initial  point  and  the 
four  corner  coordinates  of  the  box  would  be  used  to  obtain 
response  surface  values  (that  is  y  values) .  The  3 1 ze  of  the 
box  should  be  small  in  order  to  better  approximate  the 
tangent  plane  at  the  initial  point,  but  not  so  small  that  no 
effective  change  can  be  seen.  After  obtaining  the  four 
responses,  one  can  use  RSM  to  fit  an  equation  of  the  form, 
y  =  bo  +  b  1  x  1  + b2x2  +  b 1 2x 1 x 2 .  A  good  explanation  of  the  RSM 
equation  fitting  is  given  by  Meyers  (6:43-50) . 

The  following  is  how  a  2-k  factorial  design  works  on 
the  initial  point  ( x 1 0  ,  x  2  0 )  .  Let  rf  be  the  distance 
selected  in  the  xl  and  x2  direction  for  the  size  of  the  box. 


T»T 


The  four  corners  would  then  be: 

(xln.x2p.ylh)  ( x 1 p  ,  x2p  ,  yhh ) 

(  x  1  n  ,  x2n  ,  y  1  1  )  (xlp.x2n.yhl) 

where 

xln=xlO-rf 
x2n  =  x20  -r  f 
x lp  =  x  10  +  r  f 
x2p  =  x20  +  r  f 

The  first  order  equation  fitting  these  would  have 
coef  f lcients  : 


BO = ( y 1 1 + y 1 h  +  y hh  +  y h  1  )  / 4  (1) 
Bl= (yhh+yhl-ylh-yll)  4  (2) 
B  2  = (ylh+yhh-yll-yhl) /4  '3> 
B12=(yll+yhh-ylh-yhl)/4  (4) 


This  first-order  equation  approximates  the  tangent 
plane  at  yO.  B1  provides  the  gradient  in  the  xl  direction 
and  B2  the  gradient  in  the  x2  direction.  However,  if  B12  is 
not  zero,  then  there  is  interaction  and  the  surface 
approx l ma  tion  is  not  a  plane  but  a  curved  surface  Thus  »he 
B1  and  B2  slopes  could  be  misleading  if  the  interaction  is 
large  . 

The  last  exploratory  search  for  a  gradient  to  he 
examined  is  the  3-k  factorial  design.  This  design  uses 
eight  exploratory  experiments  and  fits  an  equation  to  a 
second-order  equation.  This  design  provides  more 
information  about  the  curvature  of  the  response  surface 
around  the  initial  point  The  following  is  a  description  of 
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the  3-k  design. 


The  3-k  design  uses  9  points  in  a  symmetric 


square  pattern.  In  addition  to  the  four  corner  points  of 
the  2-k  design,  the  3-k  needs  four  additional  points. 


(xln.x2p.ylh)  ( x 1 0 , x2p , y mh )  (  x  1  p  ,  x2p  ,  yhh ) 

(x In ,x20 ,y lm)  (x 10 ,x20 , ymm)  (xlp,x20,yhm) 

(  x  1  n  ,  x2n , y 1 1 )  ( x 10 , x2n , yml )  (  x  1  p  ,  x 2 n . y h 1  I 


Figure  1.  Nine  Points  of  3-K  Factorial  Design  (6:51) 

With  these  nine  points,  RSM  can  use  this  design  to  fit 
a  second-order  equation  to  this  response  surface.  The  RSM 
equation  is  of  the  form:  y  =  B0*-Blxl+Bllxl2  +  B22x2i  +  B12xlx2. 

The  coefficients  of  the  equation  are  computed  as  follows: 


BO  =  ( y 1 1 + y 1 m+ y 1 h +  y ml + ymm+ymh ♦ y h 1 +y hm+ y hh ) / 9  (6' 

B1  -  (yhl+yhm+yhh-yll-ylm-ylh)/6  !  7 ) 

B2  =  (  y  1  h  ♦  y  mh  ♦  y  hh  -  y  1  1  -  y  ml  -  y  h  1  )  /  6  (0) 

Bll  =  (yll*ylm+ylh*yhl+yhh-2»(yml+y  mm* y mh '  ' 6  ( 9 1 

B22  =  (yll=yml+yhl+ylh  +  y mh  +  yhh-2<ylm=y  mm=  y  h  m )  / 6  1 0  > 

B12=(yll=yhh-ylh-yhl)/4  ill' 


Similar  to  the  2-K  design,  B  1  provides  the  gradient  in  *  he 
xl  direction  and  B2  the  gradient  in  the  x  2  direction 

Thus,  the  literature  search  has  provided  f'Ur  ways  to 
calculate  the  gradient.  Which  of  these  is  the  most 
efficient  for  the  required  purpose ^  The  un i d i me  ns i ona 1  and 


Hooke-Jeeves  method  use  the  fewer  number  of  experiments; 
however,  both  calculate  only  one  slope  in  each  direction. 

If  there  is  the  least  amount  of  error  in  any  of  the 
responses,  it  would  affect  the  respective  slopes  greatly. 

The  2-k,  using  just  four  points,  calculates  two  slopes  in 
each  direction  and  averages  the  result  to  get  the  gradient. 
Consequently,  it  would  be  more  capable  in  dealing  with  any 
margin  of  error.  The  3-k  calculates  three  slopes  in  each 
direction  and  averages  the  result  to  get  the  best  gradient 
for  handling  noise  error,  but  it  needs  eight  additional 
points.  It  needs  four  more  points  than  the  2-k,  but  only 
averages  1  more  slope  than  the  2-k.  From  this  comparison, 
the  2-k  design  is  the  best.  It  will  be  used  in  the  program 
to  determine  the  gradient. 

Acceleration  Phase 

The  acceleration  phase  is  the  actual  climbing  up  the 
hill.  Again  it  is  desirable  to  do  this  with  as  few 
experiments  as  possible.  Many  algorithms  have  been  proposed 
on  how  to  accomplish  the  ascent  most  effectively.  Not  all 
of  the  algorithms  will  be  discussed,  merely  those  leading  up 
to  the  algorithm  used  in  the  program. 

To  begin  with,  the  initial  line  search  can  be  thought 
of  as  an  un 1 d 1  mens i ona 1  search  along  the  gradient  vector. 

The  dilemma  is  how  big  of  a  step  to  take  along  'he  vector. 
One  idea  is  to  normalize  the  slopes  to  get  a  unit  step  along 


the  vector.  One  then  takes  uniform  unit  steps  along  the 


vector.  If  the  uniform  unit  step  size  is  too  small,  it  will 
take  numerous  steps  to  reach  the  peak  of  the  vector. 
Likewise,  if  the  step  is  too  large,  the  climber  might  step 
way  beyond  the  peak.  Thus,  the  uniform  step  size  is  not  a 
very  efficient  line  search  and  ad j us  ted - s tep  line  searches 
have  been  proposed  as  an  improvement. 

Robbins- Monro  method  was  one  of  the  first  and  simplest 
improvements  over  uniform  step.  It  is  based  on  the  harmonic 
sequence  1 , 1 / 2 , 1 / 3 , 1 / 4 ,  etc.  times  the  magnitude  of  the 
dependent  variable.  The  harmonic  sequence  is  divergent  and 
the  sum  of  all  its  terms  is  infinite.  Therefore,  it 
guarantees  the  procedure  will  eventually  reach  the  peak,  no 
matter  how  far  away  it  started  (9:  162-167)  .  The  problem 
with  this  line  search  is  the  exorbitant  number  of 
experiments  necessary  to  find  the  peak,  especially  if  one 
starts  in  a  fairly  flat  region  far  from  the  optimum. 

Keston  has  devised  a  procedure  wh ich  accelerates  the 
search  more  quickly.  Instead  of  starting  with  the 
decreasing  harmonic  sequence,  Heston’s  method  starts  with  a 
uniform  step  then  shortens  the  step  size  harmonically  when 
the  peak  is  crossed  and  the  direction  of  search  reverses 
Table  I  compares  the  two  methods. 
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lerated  Peak-Seeking  Methods 


steps  12345678  total 
direction  ♦  + 

unacceler  1  1/2  1/3  1/4  1/5  1/6  1/7  1/8  1  149/280 
accelerate  1  1  1  1/2  1/2  1/3  1/4  1/5  2  17/60 


8  total 


(9 : 180) 


As  with  Robb i ns - Monr o  ,  the  Keston  method  insures  one  of 
finding  the  peak  and  it  achieves  results  more  rapidly. 
However,  it  still  takes  numerous  experiments  along  the 
vector  to  find  the  peak  of  the  vector. 

An  improvement  over  the  Keston  is  the  golden  section 
search.  It  uses  a  large  step  size  to  cross  over  the  peak. 
Once  the  peak  is  crossed  an  interval  exists  wherein  the  peak 
is  located.  The  golden  search  technique  can  now  be  used  f o 
reduce  the  interval  of  uncertainty.  Golden  search  splits 
the  interval  with  the  peak  into  two  segments  such  that  the 
ratio  of  the  whole  interval  to  the  larger  segment  is  the 
same  as  the  ratio  of  the  larger  segment  to  the  smaller. 

The  golden  search  plan  works  as  follows:  Let  the 

initial  interval  that  brackets  the  peak  be  called  d  with 
endpoints  of  zl  and  z2.  Next,  place  2  experiments  inside 
this  interval  z3  and  z4  such  that  z3=zl*0  38»d  and 
z4=zl+0.62»d.  If  the  y  response  of  z3  is  larger  than  that 
for  z4,  the  interval  of  uncertainty  is  from  zl  to  z 4 . 
Otherwise,  if  the  y-response  of  z4  is  larger  than  z  3 ,  the 
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new  interval  is  from  z3  to  z2 .  This  procedure  is  then 
continued  for  the  new  interval  of  uncertainty.  Golden 


search  will  reduce  the  initial  length  of  the  interval  of 
uncertainty  by  (O.eiS)0"1  where  n  is  the  number  of 
experiments  used.  For  example,  if  eleven  experiments  were 


used,  the  new  interval  would  be  0.008  the  size  of  the 
original  interval  (4:42-43)  .  In  addition.  Himmelblau 
recommends  a  sequential  series  of  larger  and  larger  steps 


along  the  vector  to  expedite  the  initial  bracketing  of  the 
peak.  The  golden  search  is  quite  efficient  compared  to  the 
previous  methods  and  other  interval  uncertainty  methods. 
Wilde  provides  an  excellent  comparison  of  golden  section  to 
other  interval  methods  (9:28:29) .  However,  there  is  a 


method  of  fitting  a  polynomial  to  the  points  that  is  even 
more  efficient  than  golden  search. 

The  last  un l d i mens l ona 1  line  search  that  is  examined 
and  the  one  used  in  this  computer  program,  is  the  Davies. 


Swann,  and  Campey  (DSC) -Powell  Search.  This  method  involves 
bracketing  the  peak  (DSC  portion)  and  the  fitting  of  a 
quadratic  equation  fPowell  portion)  to  interpolate  an 
estimate  of  the  peak.  G.F.  Coggins  shows  that  this 
technique  involving  the  fitting  of  a  second-order  polynomial 
through  selected  points  was  better  at  locating  the  peak  to 
within  a  specified  precision  than  the  interval  methods  such 
as  golden  section  (4:44). 
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The  Davies,  Swann,  and  Campey  (DSC)  portion  is  used  to 


bracket  the  peak.  It  involves  doubling  the  step  size  for 
each  step  until  the  peak  is  overshot.  After  the  peak  is 
overshot,  the  direction  is  reversed  and  previous  interval 
used  is  reduced  by  one-half.  This  is  used  to  obtain  one 
more  point.  This  procedure  will  give  four  equally  spaced 
points.  The  two  middle  points  y-values  are  compared  for 
optimum.  The  point  with  the  optimum  y-value  plus  the  two 
points  used  for  fitting  the  quadratic  since  the  peak  should 
be  inside  this  interval.  See  figure  2  below. 


Figure  2.  Davies,  Swann,  and  Campey  Technique 


Powell's  equation  carries  out  a  quadratic  approximation 
using  the  three  points.  The  optimum  (critical  point'  1 
found  by  taking  the  first  derivative  of  the  equation. 
Powell's  equation  (4:46)  is  as  follows: 


x-  =  -  .  5  [  (x2  J  -  x3  2  )  1 1  +  Lx3  3  -  x  1  *  >  t2  +  ( x  1  3  -  x2~ )  1 3  1 
(  x  2  -  x  3 )  t  l  +  ( x  3  - x  1  )  t  2  ♦  (  x  1  - x  2  >  1 3 
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The  equation  used  in  the  computer  program  is  of  a 
slightly  different  form.  The  derivation  of  the  computer 
equation  is  as  follows: 

Let  pi,  p2 ,  and  p3  be  the  three  points  and  tl,  t2,  and 
t3  be  their  respective  y-value.  Let  c  be  the  equal  distance 
between  the  points.  The  quadratic  equation  to  be  fitted  1 
of  the  form  y=bo+blx+b2xi  with  derivative  dy/dx=2bx+bl-Q 


implying  the  optimum  of  the  equation  is  x=-bl/2b2 


Putting 

the  three  points 

into  the 

equation 

and  solving 

simultaneously  one  gets 

bO 

+blpl+b2pl“=tl 

(  1  3  ) 

b0+blp2+b2p22=t2 

(14) 

bO 

+  blp3  +  b2p32:;t3 

i  :  5 ) 

note: 

p 1 =  p 1  p2  =  p  1  +  c 

p3  =  p  1  +  2c 

by  matrix  notation. 

bO 

b  1 

b2 

r1' 

P  1 

Pi2 

t  1  — j 

1 

p  1  ♦  c 

(  p  1  +  c  )  3 

1 2 

Li 

p  1  +  2c 

( p 1 ♦ 2c ) * 

1 3  — 1 

1-1 

P  1 

Pi2 

t  1  — | 

1  0 

c 

2cp  1  +c“ 

1 2  -  t  1  ! 

— 0 

2c 

4cpl  +  4c*“ 

1 3  -  t  1  — ^ 

i — 1 

0 

-pi  ( p  1  +  c ) 

(  c  1 1  -  p  1 

\  2  *  p  i  t  ’  '  ■  ,  — , 

!  0 

1 

2  p  1  +  c 

(  t  2  -  t  1  > 

/  c 

1 — 0 

0 

2c“ 

t  1  -  2 1 2  ♦ 

1 3  — ‘ 

note  : 

let  d=tl-2t2+t3 

and  e= -3c t 1 ♦ 4ct2 - 

c  t  3 

r—  1 

0 

0 

f  p  1  -1  *  d  + 

pl«e  +  ?cJtl  )  ''2c-  — 

0 

1 

0 

(  -  2  p  1  »  d 

» e ) .  2c- 

— 0 

0 

1 

d  /  2  c 

-- 

x  “  = 

- ! -2pl»d+e) /2c“ 

(16) 

2  (  d  /  2  c  “  ) 

x  “  - 

p  1  -  .  5  <  e  /  d  ) 

(17) 
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This  final  equation  is  used  in  the  computer  program. 

It  looks  quite  different  from  Powell’s  equation.  However, 
if  x2  and  x3  are  substituted  by  xl+c  and  xl+2c, 
respectively,  then  Powell's  equation  reduces  to  the  one 
above.  As  mentioned  earlier,  the  DSC-Powell  method  is  the 
line  search  selected  for  the  program.  This  selection  is  due 
to  its  quickness  in  finding  the  peak  and  its  accuracy. 

After  finding  this  first  point  (pi)  from  the  initial 
point  (pO) ,  one  can  repeat  the  gradient  phase  around  pi  to 
find  a  new  direction  to  proceed.  The  line  search  (DSC- 
Powell)  is  again  employed  to  find  the  next  point  (p2>  .  A 
repetition  of  gradient  and  line  searches  can  be  continued 
until  the  optimum  is  reached.  This  intuitively  attractive 
idea  of  climbing  the  steepest  path  is  known  as  the  gradient 
method,  or  the  method  of  a  steepest  ascent  (9: 107) ,  The 
advantages  of  the  steepest  ascent  method  are:  (1)  It  tends 
naturally  to  avoid  saddlepoints.and  (2)  It  will  eventually 
converge  for  any  uni  modal  function,  even  when  there  is 
appreciable  experimental  error  (9:120-121).  The  steepest 
ascent  method  is  one  of  the  two  algorithms  of  the  computer 
program.  In  the  program,  it  is  called  the  gradient/line 
me  thod . 

An  algorithm  that  accelerates  faster  than  steepest 
ascent  is  the  parallel  tangent  (PARTAN)  method.  There  are 
several  variants  of  the  PARTAN  method.  The  variant  used  is 
the  "steepest  ascent  PARTAN'.  This  method  is  also  often 


referred  to  as  'gradient  PARTAN' .  The  PARTAN  method  is  just 
like  steepest  ascent  in  finding  the  first  two  points  pO  ,  and 
pi.  After  finding  pi,  the  PARTAN  method  eliminates  the 
experimental  design  around  pi  for  the  next  direction.  The 
gradient  to  be  used  at  pi  is  the  gradient  perpendicular  to 
the  gradient  of  pO.  This  direction  will  form  a  plane  that 
is  parallel  to  the  contour  tangent  plane  of  pO  ,  where  the 
name  PARallel  TANgent  (abbreviated  PARTAN)  comes  from  The 
perpendicular  gradient  can  be  found  easily  by  swapping  the 
previous  slopes  and  reversing  the  sign  of  one  of  them.  That 
is,  if  bl  and  b2  were  the  slopes  at  pO ,  then  now  at  pi  the 
slopes  are  bl--b2  and  b2=bl.  A  line  search  is  then 
accomplished  along  this  plane  to  find  the  peak,  which  is 
point  p2.  After  finding  p2,  PARTAN  eliminates  another 
experimental  design  around  p2  for  the  next  direction. 

Instead  of  the  2-k  factorial  design,  it  connects  a  line  from 
pO  through  p2  for  a  new  gradient  direction.  A  line  search 
is  then  accomplished  starting  at  p2  and  going  along  this 
vector  direction  to  find  p3.  P3  is  the  optimum  or  very 
close  to  it.  Wh  en  the  response  surface  contours  are 
concentric  ellipsoids,  PARTAN  will  locate  the  optimum 
exactly  after  no  more  than  2k-l  unidimensional  line  searches 
(where  k  is  the  number  of  independent  variables)  (9:  1241  . 
This  means  that  point  p.3  (mentioned  above'  will  he  the  exact 
optimum  for  a  2  independent  variable  deterministic  problem. 
Thus,  after  one  initial  gradient  search  (4  experiments)  and 


three  line  searches  (about  5  experiments  each)  ,  PARTAN 


* 


* 


locates  the  optimum.  Whereas  the  grad  1 ent / 1 1 ne  uses  two 
gradient  searches  (4  experiments  each  and  two  line  searches 
(about  5  experiments  each)  per  zigzag.  Therefore,  one  can 
conclude  the  PARTAN  algorithm  is  indeed  a  more  efficient 
method  for  certain  quadratic  problems. 

Even  when  the  contours  are  not  precisely  elliptical, 
PARTAN  has  certain  ridge  following  properties  which  make  it 
attractive  especially  when  the  ridges  are  straight  (10:323) . 
In  addition,  PARTAN  will  work  perfectly  in  two  dimensions 
for  any  radially  similar  contours  since  the  property  of 
parallel  tangents  works  for  these  (9:144).  Even  for  other 
non - e 1 1 l pso l da  1  surfaces,  PARTAN  can  still  work.  It  is  just 
that  PARTAN  will  generally  not  be  right  at  the  optimum  after 
one  cycle,  but  this  does  not  prevent  starting  over  again 
using  point  p3  as  the  beginning  of  another  PARTAN  search. 

The  geometric  reason  PARTAN  works  (finds  the  top  of  the 
hill  with  so  few  line  searches)  can  be  simply  explained 
using  a  contour  plot  of  the  response  surface.  See  figure  3. 
Let  a  point  pO  be  randomly  selected  and  a  line  be  dra wn  from 
pO  to  the  center  of  the  ellipse,  p3.  One  will  notice  that 
this  line  p0p3  intersects  each  contour  ring  at  the  same 
angle.  Next,  the  contour  tangent  planes,  t ( l ) s ,  are  drawn 
at  the  point  of  each  of  these  intersections.  One  will 
observe  the  planes  are  all  parallel.  Also,  the  point  of 
Intersection  with  the  contour  ring  is  the  optimum  point 
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along  the  plane  line  for  each  tangent  plane.  Next,  the 
gradient  vector  g  at  pO  is  perpendicular  to  the  contour 
tangent  plane,  tO  ,  at  pO  .  Thus,  perpendicular  to  all  the 
other  contour  tangent  planes  drawn. 

The  PARTAN  method  described  above  would  follow  the 
darkened  path  in  figure  3.  Starting  at  pO,  PARTAN 
calculates  the  gradient  vector.g.  It  goes  along  this  vector 
to  the  vector  peak,  pi  It  then  mo  ves  on  a  vector 

perpendicular  to  g  at  pi  to  the  vector  peak,  d  2  It 
connects  pO  to  p  2  and  follows  this  vector  to  its  peak,  p  3  . 
This  described  method  will  be  called  the  PARTAN'line  method 
for  the  rest  of  this  paper  and  in  the  computer  program 
There  is  still  another  improvement  to  the  search 
method.  Faix  proposes  an  efficient  improvement  upon  the 
PART  AN /line  method.  Faix  states. 

The  method  only  works  exactly  for  perfect 
quadratic  response  surfaces  with  no  noise.  However,  it 
will  be  shown  to  be  relatively  robust  against  many 
types  of  imperfection,  and  thus  a  good  methodology 
choice  [1:180]. 

This  improvement,  to  be  called  the  PARTAN/FAIX  method, 
eliminates  the  line  search  from  p2  to  p3.  The  PARTAN 'FA  IX 
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me thod  calculates  the  distance  from  p  2  to  p  3  by  using  the 


results  from  points  pO ,  pi,  and  p2  to  find  the  eccentricity 
of  the  ellipse.  The  eccentricity,  e,  measures  the  stretch 
of  the  ellipse.  The  eccentricity  of  an  ellipse  is  equal  to 
c/a  in  figure  4  and  ranges  from  0  to  almost  1. 


Figure  4.  Eccentricity  of  an  Ellipse 

The  eccentricity  is  zero  for  a  circle  and  approaches  one  as 
the  major  diameter  increases  in  ratio  to  the  minor  diameter 
(10:379-399).  For  a  quadratic  equation, 
y  =  bo  +  b  1  x  1  + b2x2  +  b  1  1  x  1  2  +  b2 2 x 2 ~  ,  one  can  relate  the 

coefficients,  bll  and  b22,  to  e.  The  square  root  of 
(b22/bll)  is  equal  to  the  a/b  in  figure  4.  Therefore,  if 
b22  is  greater  than  bll,  then  e  is  equal  to  the  square  root 
of  (l-bll/b22)  and  if  bll  is  greater  than  b22,  then  e  is 
equal  to  the  square  root  of  (b22/bll-l).  In  agree  me  n  t  with 
Faix's  notation,  b  2  2 / b  1  1  will  be  called  the  variable  c 
(1:186). 
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The  variable  c  can  be  found  geometrically  using  the 


points  pO,  pi,  and  p2  of  a  PART AN  search.  Using  figure  3, 

c  =  (  1  +  mo  *  r  )  /  (  mo  »  r  -  (  mo  )  *  *  2  )  (18) 

where 


mo  = 

r4/r5,  the  known  slope 

of  line  r3  bet  we  e  n 

pO  an  p3 

r  = 

rl/r2,  the  ratio  of 

lines  r 

1  and  r2 

r  1  = 

the 

distance 

be  tween 

pO 

and 

P  1 

r  2  = 

the 

distance 

bet  we  e  n 

P  1 

and 

P  2 

r  3  - 

the 

d 1 s  tance 

be  tween 

pO 

and 

P  2 

r  4  = 

the 

distance 

x 22  - x20 

r  5  = 

the 

distance 

x  12-x 10 

(  1  :  182 ) 

Using  the  variables  c,  mo ,  and  r  3  .  Faix  derives  the 
length  of  the  acceleration  step  between  p2  and  p3  in 
multiples  of  r3 .  The  length  between  p2  and  p3  equal  r3*assu 
where  (1:182) 

assu= [ (mo) *« (c- 1 ) *«c ]/ [ I + (c*mo) 2 ] J  (19) 

There  is  an  equally  efficient  method  to  the  PART  AN /FAIX 
method.  One  may  notice  that  instead  of  using  the  parallel 
plane,  t3 ,  any  of  the  other  parallel  planes  would  have 
worked  for  PARTAN.  Thus,  instead  of  doing  a  line  search 
from  pO  to  find  pi,  choose  any  distance  to  place  pi  from  pO . 
A  unit  step  of  1  is  offered  as  an  option  in  the  computer 
program.  Then  do  a  line  search  perpendicular  to  find  p2  . 
Connect  pO  and  p2  and  do  a  line  search  in  this  direction  to 
find  p3.  This  would  entail  just  two  line  searches 
comparable  to  the  PARTAN/FAIX  method. 


s  V 


Second-Order  Exploratory  Phase 


Both  the  beginning  and  the  end  of  a  search  involve 
local  exploration  The  beginning  being  a  simple  linear 
study  near  an  arbitrary  point  to  get  to  the  gradient 
direction.  At  the  end  of  the  search,  a  nonlinear 
exploration  in  the  vicinity  of  the  optimum  is  accomplished 
to  insure  the  optimum  was  found  (9:75) .  This  end 
exploration  may  actually  find  a  point  nearby  that  is  better 
then  the  optimum  found  by  the  algorithm  This  could  be 
caused  by  error  in  the  simulation  or  testing  process  In 
addition,  this  final  exploratory  phase  will  show  the 
behavior  of  the  response  surface  near  the  optimum. 

A  3-k  factorial  design  is  used  to  provide  the  second  - 
order  equation  fit.  This  design  was  discussed  under  the 
gradient  phase.  The  difference  now  is  once  the  coefficients 
are  found,  they  are  fitted  into  a  derived  equation  for  the 
critical  point. 

x  1  f  =  x  1  0 ♦ [  ( - b2  *  b 1 2 )  +  (  2  »  b  1  l«bl)  ]/[  (bl2»b!2)-(4*bll*b22)  ]  (20) 
x2f-x20+[  (  -  b  1  *  b  1  2  )  ♦  ( 2  *  b 1 1  *  b2 )  ]/[(bl2»bl2)-(4*bll»b22)  ]  (21) 

If  (4*bli*b22)  is  less  than  (fcl2*bl2)  ,  then  this 
point  is  a  saddle  point.  Otherwise,  if  bll  is  less  than  0 
and  b22  is  less  than  0,  then  the  point  is  a  maximum.  If  bll 
and  b22  are  positive,  then  the  point  is  minimum.  Thus,  the 
coefficients  bll  and  b22  describe  the  shape  of  the  response 
surface  around  the  optimum. 
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Summary 

This  chapter  reviewed  the  literature  and  methodology 
that  lead  up  to  the  writing  of  the  compute”  program.  The 
next  chapter  describes  the  actual  computer  program  and  how 
1  t  wo  r  k  s 
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III.  Program  Description 


This  chapter  describes  the  flow  of  the  computer 
program  The  complete  program  is  contained  in  Appendix  A. 
The  program  is  written  in  FORTRAN  on  a  VAX  computer  but 
could  be  transferred  to  a  microcomputer  for  further  use 

Before  running  the  program,  the  subroutine  SIM  must  be 
modified  in  t  wo  wa  y s 

1)  Line  18  of  the  subroutine  must  reflect  whether  the 
problem  is  a  minimization  or  ma  ximization.  This  is  done  tv 
removing  the  letter  c  in  the  first  colu  mn  of  line  18  for  a 
minimization  and  ensuring  the  letter  c  is  in  place  for  a 
maximization  The  letter  r  comments  out  line  18  for  a 
maximisation  problem 

<  7  >  The  problem  simulation  must  be  loaded  into  the  SIM 
subroutine  starting  at  line  14.  If  the  simulation  or 
experimentation  is  to  be  run  externally  of  the  program, 
column  1  in  line  14  gets  a  Letter  c  added  and  column  1  in 
lines  1 c'  .  16,  and  1  ”  get  the  letter  c  re  mo  ved  . 


Main  Program 

The  main  program  is  called  SEARCH  A  f 1  ' w  1  i a  g  r  a  m 


SEARCH  is  sho  wn  in  figure 


The  program  has  rumen 


interactive  options  for  the  user  *  a  1  .  - w  t  he  user  r  he 


freedom  to  wo  rk  a  variety  of  pr't.e  ms  However  •'he  program 


is  mainly  designed  to  efficiently  optimize  quadratic 
problems 

SEARCH  consists  of  repetitions  of  asking  the  user  to 
make  a  choice  and  then  calling  the  subsequent  subroutines 
and  showing  the  outcome  for  that  choice.  This  allows  the 
program  to  step  along  from  point  to  point  to wa rd  the 
o  p  1 1  mum , 

The  program  begins  by  asking  the  user  for  the  point  0 
xl  and  x2  starting  coordinates.  It  then  calls  the 
subroutine  SIM,  which  gives  the  y  response  for  these  inputs 
and  the  main  program  writes  these  values  to  the  screen  and 
output  file.  At  this  point,  it  automatically  calls  the 
subroutine  TWOK.  TWOK  does  a  2-k  factorial  design  and  RSM 
fit  to  find  the  gradient  directions.  If  the  linear  equation 
is  of  a  flat  surface,  thus  having  no  gradients,  the  mam 
program  will  end  for  there  is  no  direction  to  climb  at  this 
point.  Otherwise,  the  ma in  program  will  write  to  the  screen 
and  the  output  file,  the  nor  ma lized  gradient  directions. 

This  gradient  direction  is  the  best  direction  for  climbing 

Next  the  main  program  will  prompt  the  user  for  how  far 
to  travel  in  this  gradient  direction.  The  two  options  are: 
one  unit  step  or  to  the  peak  in  that  direction.  The  first 
option  should  be  used  only  if  the  response  surface  is 
ellipsoidal.  The  second  choice  might  be  used  with  PARTAN  m 
the  steepest  ascent  method.  With  the  choice  made,  the 


program  will  calculate  point  1  and  write  the  location  of 
point  1  to  the  screen  and  output  file. 

Proceeding  from  point  1,  the  main  program  calculates 
the  gradient  perpendicular  to  the  last  gradient  direction 
and  writes  it  to  the  screen  and  output  file.  It  then  calls 
subroutine  LINE  to  find  the  peak  along  that  gradient  line. 
After  calling  SIM,  it  wr ites  the  calculated  point  2  to  the 
screen  and  output  file. 

At  point  2,  the  user  decides  to  use  the  PARTAN  method 
or  continue  the  grad i ent/ 1 i ne  method.  If  option  1  (gradient 
/line)  is  selected,  then  the  perpendicular  gradient  is 
calculated  and  written.  The  program  calls  subroutine  LINE 
to  find  the  peak  in  this  direction.  This  is  followed  by  the 
subroutine  SIM.  The  grad/line  point  3  is  then  written  to 
the  screen  and  output  file.  If  option  2  (PARTAN)  is 
selected,  then  the  subroutine  PARTAN  is  called  to  calculate 
the  gradients.  These  two  gradients  are  written  to  the 
screen.  With  the  PARTAN  directions,  the  program  offers  the 
user  the  option  of  doing  a  line  search  or  a  FAIX  jump  to  'he 
next  point.  If  a  one  unit  step  was  selected  at  point  0, 
then  a  line  search  must  be  selected.  Otherwise,  the  second 
option  (FAIX  method)  is  the  most  efficient  and,  if  selected, 
the  main  program  calls  the  subroutine  FAIX  to  compute  the 
PARTAN/ FAIX  point  3. 

Next,  the  main  program  offers  the  user  to  choose  which 
of  the  previous  three  options  he  wants  to  use  for  point  3. 


This  is  included  in  case  more  than  one  option  was  selected. 
Point  3  is  then  written  to  the  screen  and  the  output  file. 

Finally,  the  mam  program  asks  the  user  if  he  wants  to 
exit  at  this  point,  repeat  the  whole  process  again  using 
point  3  as  the  new  initial  point,  or  to  do  a  3-k  factorial 
design  and  RSM  to  better  locate  the  final  point.  If  the 
user  is  confident  the  surface  is  ellipsoidal  and  there  was 
little  error  in  the  input  values,  then  one  should  be  at  the 
optimum  and  exiting  is  the  correct  choice.  If  the  user  is 
sure  the  optimum  has  not  yet  been  reached,  maybe  due  to  the 
complexity  of  the  surface,  then  repeating  the  process  again 
would  produce  the  better  answer.  If  the  user  feels  close  t 
the  optimum,  but  point  3  is  slightly  off,  then  3-k  design 
with  RSM  will  help  to  find  the  final  optimum 


T MLK._S.uh  ra  u  t_i.n  e 


£ 


« 


This  subroutine  fits  a  2-k  factorial  design  around  a 
point  and  uses  RSM  to  fix  an  equation  to  the  four  points. 

The  coefficients  of  the  equation  are  used  for  the  gradients. 
Figure  6  is  the  flow  diagram  for  the  TWOK  subroutine.  The 
subroutine  begins  by  stating  the  radius  of  the  2-k  factorial 
design  and  asking  the  user  if  he  would  like  to  change  the 
radius  size.  After  the  radius  size  is  determined,  the 
coordinates  of  the  four  corner  points  are  calculated.  With 
these  coordinates,  the  subroutine  SIM  is  called  and  the 
response  values  found.  Next,  the  maximum  y-value  of  the 
four  points  is  called  the  variable,  m.  This  variable  is 
compared  to  the  initial  point  response,  yO .  If  vO  is  larger 
than  m,  then  there  appears  to  be  no  uphill  direction  from 
the  initial  point.  Therefore,  the  subroutine  RSM  is  called 
for  an  exploratory  search  of  the  optimum  within  this  area 
Otherwise,  the  four  y-responses  are  used  to  calculate  the 
coefficients  of  the  first-order  equation.  Of  these  four 
coefficients  bO,  b 1  ,  b  2 ,  and  bl2,  bl  and  b2  are  used  as  the 
xl  slope  and  x2  slope,  respectively.  If  both  of  these 
values  are  zero,  then  there  is  no  slope  in  this  area  and  the 
program  will  print  'Try  a  new  starting  point.'  and  end 
Before  bl  and  b2  are  passed  back  to  the  main  program,  the 
subroutine  normalizes  their  value.  The  control  then  returns 
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to  the  main  program. 
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Figure  6.  Flow  Diagram  of  TWOK  Subroutine 


LINE  Subroutine 

This  subroutine  finds  the  peak  in  a  vector  direction 
using  the  DSC-Powell  algorithm.  Figure  7  is  the  flow 
diagram  for  the  LINE  subroutine.  The  subroutine  first 
calculates  a  point  that  is  2  units  in  the  gradient  direction 
from  the  starting  point.  After  calling  SIM  to  get  the  y - 
response,  it  checks  to  see  if  the  response  wa  s  an  increase 


III  VS 


over  the  initial  response.  If  it  was  not  an  increase,  the 
subroutine  reverses  the  gradient  direction  and  calculates  a 
point  2  units  in  the  other  direction.  It  again  calls  SIM  to 
get  the  y-response.  If  this  too  was  not  an  increase,  the 
program  curve  fits  these  three  points  to  find  an  optimum. 

If  either  direction  had  been  an  increase  response,  the 
subroutine  wo uld  double  the  step  size  and  calculate  the  next 
point.  It  would  continue  doubling  the  step  size  until  it 
has  either  gone  10  steps  ( to  prevent  a  continuous  loop)  or 
got  a  response  that  wa s  a  decrease  from  the  previous  step. 
Once  it  gets  a  y-response  that  is  a  decrease,  it  cuts  the 
step  size  in  half  and  reverses  the  vector  direction.  This 
gives  four  equally  spaced  points.  The  subroutine  compares 
the  2  middle  responses.  It  uses  the  point  with  the  larger 
response  and  the  points  on  both  sides  of  it  to  curve  fit  an 
equation.  The  first  derivative  of  this  equation  is  used  to 
find  the  optimum  point  along  the  vector.  This  optimum  point 
is  then  returned  to  the  main  program. 


by  asking  the  user  how  many  repetitions  of  the  problem  is  to 


be  accomplished  at  these  values.  The  default  is  1 
simulation.  It  then  loops  through  the  simulation  or  user 
input  the  required  number  of  repetitions.  The  subroutine 
then  averages  the  responses  to  get  one  value  to  pass  back  to 
the  program.  Also,  by  enabling  line  18  the  program  can  run 
minimization  problems  by  doing  a  negative  maximization. 


Figure  8.  Flow  Diagram  of  SIM  Subroutine 


PABTAN.  Sub.rautin.fi 

This  short  subroutine  takes  the  coordinates  of  point  0 
and  point  2  and  calculates  the  slope  between  these  two 
points.  This  slope  is  then  normalized  and  passed  to  the 
main  program  to  be  used  as  the  next  gradient  direction. 
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FA IX  Subroutine 

This  subroutine  takes  the  coordinates  of  points  0,  1, 

and  2  and  calculates  the  distances  between  these  points.  It 
also  calculates  the  slope  between  point  0  and  point  2.  It 
then  computes  the  ratio  of  the  distance  between  pO  and  pi 
over  the  distance  between  pi  and  p2 .  With  these  values,  the 
subroutine  determines  values  for  the  variables  c  and  assu. 
Finally,  it  uses  assu  and  the  distance  between  pO  and  p2  to 
estimate  the  location  of  the  p3  coordinates.  After  calling 
SIM  subroutine  to  get  the  y-response  for  point  3,  it  returns 
to  the  main  program. 

RSM_.Su.br  out  me 

This  subroutine  fits  a  3-k  factorial  design  around  a 
point  and  uses  RSM  to  fit  a  second-order  equation  to  the 
nine  points.  The  first  derivative  of  the  fitted  equation  is 
used  to  find  the  critical  point  in  this  area.  The  second 
derivative  test  is  then  used  to  determine  whether  the 
critical  point  is  a  maximum,  minimum,  or  a  saddle  point. 

Summary 

This  chapter  described  the  procedures  of  the  SEARCH 
program.  It  looked  at  the  contents  of  each  subroutine 
and  the  flow  of  the  ma in  program.  The  next  chapter  discusses 
how  we  11  the  program  wo r k s . 
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IV.  Validation.,  of  Program 


A  validation  of  the  program  will  be  accomplished  by 
taking  a  known  quadratic  equation  and  comparing  the  SEARCH 
program  output  with  the  known  mathematical  values.  See 
figure  9  for  the  output  file.  The  equation  is  y  =  10- 

5  (  x  1  ♦  1  )  2  -  1  5  (  x 2  -  2 )  3  An  initial  point  of  (6,9)  will  be  used 

to  start  the  program. 

The  first  check  is  of  the  gradient  around  the  initial 
point.  The  SEARCH  program  obtains  gradients  of  (- 
0.3162278,-0.9486833).  The  first  derivatives  of  the 
equation  are  (-10(x+l)  ,-30(x-2))  and  at  the  initial  point 
(6,9)would  yield  (-70,-210)  as  the  gradients.  If  (-70,-210) 
is  normalized,  one  gets  exactly  the  same  values  as  the 
SEARCH  program  obtains.  Thus,  the  TW0K  subroutine  does 
obtain  accurate  gradients. 

The  next  check  is  of  the  line  search  to  point  1.  The 
SEARCH  program  used  four  steps  to  get  to  point  1.  It  used 
three  steps  before  it  passed  the  vector  peak  and  one  reverse 
step  to  evenly  space  the  points.  After  curve  fitting,  the 
program  obtained  point  1  as  <3.5,15).  This  point  1  is  7  q 
unit  gradient  steps  from  point  0  with  a  y  -  value  of  -  9  5  .  0  0  3  . 

To  check  the  accuracy  of  this  line  search,  one  can  ''heck  .  8 
and  8.0  unit  gradient  steps.  This  gives  v -  values  of  -95.156 
and  -95.1245.  respectively  Thus,  the  line  search  was  quite 
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accurate  in  selecting  the  peak  value  along  the  gradient 
vector . 

After  point  1 ,  the  program  computes  gradients  that  are 
perpendicular  to  the  first  set  of  gradients.  One  can  see  by 
inspection  that  the  new  gradients  (-0.949,0.316)  are 
perpendicular.  Thus  the  calculation  is  correct. 

Another  line  search  is  accomplished  after  just  three 
steps  to  find  point  2.  The  SEARCH  program  after  three  steps 
finds  point  2  as  (-0.25,2.75)  with  a  y-value  of  -1.250020. 
Using  3.9  and  4.0  gradient  steps  to  check  the  accuracy,  one 
gets  y-value s  of  -1.2657  and  -1.26336,  respectively.  Thus 
the  line  search  is  very  accurate  again. 

Next ,  the  program  calculates  the  PARTAN  gradient  as  <- 
0.707,-707).  This  is  the  slope  between  point  2  and  point  0 
(  2 . 75 -9  .  -  .  2 5 -6 )  and  are  the  same  values  once  normalized 

Finally,  the  FAIX  subroutine  calculates  point  3  as  - 
1,2)  with  a  y-response  of  10.  This  can  be  checked  by  the 
format  of  the  equation  as  the  actual  optimum.  The  RSM 
subroutine  is  also  ran,  but  shows  it  is  unable  to  improve 
the  o  p  t l mum . 

This  demonstration  by  example  has  shown  the  accuracy  of 
the  subroutines  that  make  up  the  program.  The  output  from 
other  problems  are  contained  in  Appendix  B. 

The  efficiency  of  the  SEARCH  program  is  evident  bv  ' he 
few  simulation  runs  required  The  above  problem  needed  a 
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Problem: 


5  (  x  1  ♦  1  )  *  *  2  -  15 ( x  2  -  2  )  »  »  2 


point  0  =  6 

5 . 500000 

5 . 500000 

6 . 500000 
6 . 500000 

bO  .  b  1  ,  b2  ,  b  1  2 
-975 . 0000 
slope  in  xl 
slope  in  x2 
5 . 367545 
4  .  102633 
1.572811 
2 . 837723 
point  1  = 


6 . 000000 

8 . 500000 

9 . 500000 


8 . 500000 

9 . 500000 

are 

-35 . 00000 
direction  is 
direction  is 
7 . 102633 
3 . 307900 
-4 . 281567 
-0 . 4868333 
3 . 500000 
e  for  xl  is  < 
e  for  x2  is  -( 


the  new  slope  lor  xl  is  u 
the  new  slope  for  x2  is  -0 
5.397367  0.8675440 

1  . 602634  2  .  1  32455 

-2 .  192  100  3 . 397366 

-0.2947329  2.764911 

point  2  is  -0.2500001  I 

partan  slopes  are  -0.7071068 
r 1 . r2 , r3 , r , mo  are 

7.905694  3.952848  8. 

c=  3.000001 
assu:  0.1 200000 

faix  point  3  is  -1.000000 
point  3  is  -1.000000 


9.000000 

-835 . 0000 

-  1045 . 000 
-905 . 0000 
-1115. 000 

-  105 . 0000 
-0.3 162278 
-0 . 9486833 

-  583 .2812 
-145.8434 
-614. 9680 
-156.4057 

1 . 500000 
0 . 9486833 
0 . 3 162278 


0000 


0 . 00000 COE+OO 


00003 


-2  13. 8684 
-24  .  13  168 
-26 . 39499 
-  1  .  263333 
750000 


-  1  .  249993 


8 . 838835 


-0 . 707  1068 


2 . 000000 


1  . 000000 


-  1  . 500000  1  . 500000 

-  1  . 500000  2 . 000000 

-  1  . 500000  2 . 500000 

-  .  ^00000  1  .  500000 

-  1  . 000000  2 . 000000 

-  1  . 000000  2 . 500000 

-0.5000001  1.500000 

-0.5000001  2.000000 

-05000001  2.500000 

bO.bl  , b  2 , b 1  1  ,b22,bl2  = 

6.666667  1  .  1  1 2 6 2 0 0 E - 06  3.655' 

1 . 1 920929E-07 

the  final  point  is  -1.000000 
(x 1 f , x2f )  is  a  maximum  point. 


2 . 000000 
2 . 000000 

4 . 999996 
8 . 749999 
5 . 000003 
6 . 249997 
1 0 . 00000 
6 . 250004 
4 . 999998 
8 . 75000  1 
5 . 000005 

65575  1 5E  -  06 


00000 
10 . 00000 


250000 


2  .  000000 


750000 


0  0  0  0  0 
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mere  twelve  simulations  to  find  the  optimum  plus  eight 
additional  to  run  the  RSM  accuracy  check  at  the  end. 

A  comparison  of  this  program  to  other  similar  programs 
is  not  realistic  Most  other  programs  are  written  to  solve 
complex,  nonlinear  equations.  Reklaitis  compares  several  of 
these  methods  and  algorithms,  and  was  unable  to  determine  a 
superior  me  thod  (7:60,120). 

Summary 

This  chapter  used  an  example  problem  as  a  validation  if 
the  program.  It  showed  the  accuracy  of  the  subroutines  and 
efficiency  of  the  program.  The  next  chapter  recommends 
further  enhancements. 


Conclusion  and  Re commend at  ions 


V  . 

The  overall  objective  of  this  research  was  to  develop 
an  interactive,  user-friendly  computer  package  that  allows 
one  to  locate  the  optimum  response  of  an  unknown  objective 
function  in  a  minimum  number  of  experimental  trials 

Subobiectives  were ' 

'1)  Show  that  the  techniques  chosen  were  the  most 
efficient  using  the  least  number  of  trials  as  the  measure  of 
ef  feet  ive ness 

(2)  Verify  that  the  program  can  find  the  optimal 
response  to  a  sample  problem. 

This  research  effort  accomplished  these  objectives. 

The  program  is  user-friendly  and  solves  the 
o p 1 1  mi za 1 1 onprob 1 em  in  a  very  efficient  number  of  trials 
It.  in  addition,  provides  the  flexibility  to  solve  even  more 
complex  problems  than  just  quadratic  surfaces 

Recommendations  for  further  enhancements  would  be  to 
expand  the  number  of  independent  variables  that  the  program 
can  handle.  This  would  increase  the  base  of  problems  the 
program  can  solve.  Another  enhancement  would  be  ♦ -  enable 
the  program  to  Incorporate  constraints  This  wv.l  i  f  r  :  a  1  e  n 
its  adaption  to  real  wo  rid  problems  Finally.  '•  r,  e  : 
user- friendly  enhancement  would  be  to  add  ir,  s  own  Graphic 
display  of  the  response  surface 
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program  search 

integer  g,i,j,k,l,m 

real  x 1 0 , x20  ,  yO  ,  b  1  ,  b2  ,  r  1 

real  xll ,x21 ,yl ,xl2,x22,y2,xl3,x23.y3 

real  1'  1  3  ,  f  23  ,  f  3  ,  1  1  3  ,  1  23  ,  1  3 

real  x 1 f , x2f . y f 

open  17 .file-' output '  , s  tatus  = ’ unknown ’  ) 

c - this  portion  interactively  gets  the  starting  point. 

print  »,  ’Input  the  starting  point.  This  should' 
print  *  ,  'be  your  best  guess  of  the  optimum.' 
print  *,  'Enter  the  xl  coordinate, 
read  * , x 1 0 

print  *  ,  'Enter  the  x2  coordinate.’ 
read  *  .  x  2  0 

call  SIM  (xl0,x20,y0) 

print  * ,  'The  value  for  yO  ls’.yO.’.' 

2  0 0  wr l te  17 , * )  point  0=’,xl0.x20.y0 


110  print  ».  'Around  this  initial  point  ixlO  ,x20)  a  2K 
factorial ' 

print  ».  'design  is  accomplished  to  get  the  ' 
print  «,  'gradient  (direction  of  ascent*.' 

c - this  subroutine  calculates  the  first  gradient 

call  TWOK  < x 10 , x20 , yO , b  1  . b2 i 

c - this  if  statement  is  for  flat  surface 

if  C(bl.eq.0).and.!b2.eq.0))  then 
go  to  1000 
end  if 

print  *.  ' 3  lope  in  xl  direction  is'.bi 

print  *.  ’slope  in  x  2  direction  l s  '  .  b  2 

wr i te ! 7 , « i , ' slope  in  xl  direction  is'.bi 
wr 1 te (7 , • )  slope  in  x  2  direction  l s ’  .  b  2 

print  » .  'Enter  1 )  to  go  1  unit  step' 

print  »,  'Enter  2  >  to  do  a  line  search’ 

read  *  ,  j 

if  ( j  . eq .  I )  then 
x  1  1  =  x  1  0  ♦  b  1 
x  2  1  r  x  2  0  +  b  2 

call  SIM  ( x i  1  . x  2  1  .  y  1  1 
go  to  300 
end  if 

c -  --this  subroutine  'ires  a  line  search  tor  the  next 

point  . 
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print  #,  Using  the  gradients,  a  line  search  will  be 
print  ».  ’accomplished  to  find  the  peak  in  this 
direction ' 

call  LINE  (xl0,x20,y0,bltb2,xll,x21,yl,i) 

c - this  if  statement  is  for  lines  that  never  peak 

if  (i.eq.10)  then 
go  to  1000 
end  1  f 

300  print  *,  'The  xll  coordinate  is’ ,xll 
print  »,  'The  x21  coordinate  l s' ,x21 
print  *  ,  'The  value  for  yl  is '  ,yl  ,  '  .  ' 
wr  ite(7,»)  ,  'point  1  =’txll,x21,yl 

print  *,  'Do  you  want  to  do  a  2-k  design  to 
print  »,  'get  the  next  gradient  or  use  the 
perpendicular ' 

print  »,  'to  the  line  search, 
print  »,  ' 1 )  perpendicular’ 

print  »,  ’2)  2-k  factorial  design' 

read  * , g 

1 f ( g  . eq .  1 ) then 

r  1  =b  1 
b  1  =  -b2 
b  2  =  r  1 
go  to  25 
end  if 


■this  step  gets  a  gradient  at  this  location, 
print  *,  'Another  2-k  factorial  will  be  done 
print  »,  'to  find  gradients  from  this  point 


cal  1  TW0K  ( x 1  1  , x2 1 
if  (  (  b  1  eq.0). and 

go  to  1000 
end  1  f 

print  »  .  'slope 
print  *  ,  'si  ope 
wr  l  t  e  (  7  ,  • )  ,'the 
wr  ite  (7  ,  *)  ,  'the 


y  1  ,  b  1  .  b  2  ) 
(  b'2  .  eq  .  0  ) 


then 


i  n 
l  n 


xl  direction  is',bl 
x2  direction  l  s  '  .  b  2 
new  slope  for  xl  is’ , b 1 
new  slope  for  x 2  is',b2 


: 


•this  step  does  a  line  search  for  the 
print  *  ,  'A  line  search  will  be  done 
print  *,  'gradients  from  point  1.' 
call  LINE  (xl 1  ,x21  ,yl  ,bl  .b2.xl2.x22 ,y2 
if  (  l  .  e  q  .  1  0 )  then 


next  point, 
using  these 

l  ) 


go  to 
end  if 
print  * 
print  » 
print  * 
wr  i  t,e  (7 


1000 

'The  x  1  2 
'The  x  2  2 
'The  value 
» l , ' point  2 


coordinate  is’ ,xl2 
rnnrd  i  nat.p  is'  x2? 
for  y2  i s '  ,  y  2  , 
is'  , x 1 2 , x  2  2  ,  v  2 
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c - the  next  step  lets  one  choose  between  another 

grad  lent/ line 

c  search  or  the  shorter  partan . 

print  »,  'Do  you  want  to  do  another  gradient/line 
search  or ’ 

print  *,  'the  partan  method  or  both7’ 

print  *  ,  'The  partan  connects  point  0  and  2  for  the 
g rad  lent.  ' 

print  *,  '  1)  yes  -  g r ad l en t / 1 l ne ’ 

print  «,  '  2)  yes  -  partan’ 

print  »,  ’  3)  ves  -  compare  both' 

print  *,  'enter  1  , 2  .  or 3  .  ’ 
read  *  ,  k 

if  (k.ne.2)  th<-r 

call  TWOK  1 x 1  , x22  ,  y2  .  b  1  ,  b2 ) 

print  *  ,  '  1  >io  x  b  1  =  ’  ,  b  1 
print  *  ,  ’ t  wo  k  b  2  = ’  , b  2 
write(7,«) , ' twok  slopes  are' ,bl ,b2 

call  LINE  ( x 1 2 , x22  ,  y2  ,  b  1  ,  b2  ,  x  1  3  ,  x23  ,  y3  ,  l  ! 
print  *,  ’ grad l ent/ 1  l  ne  new  point  location  is' 

print  *  ,  x 1 3 , x23  ,  y3 

write  (7  ,  *  )  ,  '  grad  /  1  me  point  3  is’.xl3,x23,y3 
end  if 

if  (k.ne.l)  then 

call  PARTAN  ( x  1  0  .  x20  ,  x  1  2  .  x22  .  b  1  ,  b2 ) 
print  *,  'partan  bl  =’,bl 

print  »,  'partan  b2  =’  ,b2 

write(7,»)  ,  'partan  slopes  are'  ,bi ,b2 

print  *,  'do  you  want  to  do  a  line  search  or  the 
faix  method7' 


print 

# 

'  1  ) 

yes  -  line  search’ 

print 

# 

'  2 ) 

yes  -  faix  method’ 

print 

# 

’  3  ! 

yes  -  compare  both 

print 
read  » 

# 

,  1 

'enter  1  , 2 , or3  .  ’ 

if  ( 1 . ne . 1 ) then 
call  FAIX 

(xi0.x20.xll  ,  x  2 1  ,  x  1  2  ,  x  2  2  .  b  1  ,b2.fl3.f23,f3) 

print  *,  'the  FAIX  method  calculated' 
print  »  ,  f  1  3  ,  f  23 , f  3 

wr  l  te (7  ,  *)  ,  '  faix  point  3  is’,fl3,f23,f3 
end  l  f 

if  (  1  .  n  e  .  2 )  then 

call  LINE  (xl2.x22.y2.bl.b2, 113, 12 3, 13,i) 
print  »,  'the  line  search  found' 
print  »,  113,123,13 

write(7,*)  ,  'partan/line  point  3  is'.  113.  12 3.  13 
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end  1  f 
end  if 

print  *,  'which  do  you  want  to  use9' 

print  »,  'enter  1  to  use  partan/FAIX  values  (if 

used )  ' 

print  *  .  f 1 3  ,  f  23  ,  f  3 

print  *  ,  'or  2  for  partan/1  me  values  (if  used) 
print  * ,  113,123,13 

print  »,  'or  3  to  use  g r ad i en t / 1 i ne  values  (if 

used)  ’ 

print  » ,  x  1  3  ,  x  2  3  ,  y  3 

read  *  ,  m 
if  (m.eq.l)  then 
x  1  3  =  f  1  3 
x  2  3  =  f  2  3 
y  3  =  f  3 

else  if  (m.eq.2)  then 
x  1  3  =  1  13 
x  2  3  =  12  3 
y3  =13 
end  if 

print  *,  'The  xl3  coordinate  is’  ,x!3 
print  *,  'The  x23  coordinate  is’.x23 
print  *  ,  ’The  value  for  y3  l s ’  ,  y 3  ,  ’  .  ’ 
write(7,*) ,  ’point  3  is’,xl3,x23,y3 
print  »,  ’which  do  you  want  to  do9’ 
print  *,  ’1)  quit/exit' 

print  *,  ’2)  repeat  process  using  point  3  as  lnitia 

point  ’ 

print  »,  ’3)  3-k  factorial  design  and  RSM’ 

read  * , h 

if  (h  .eq.  1  ) then 
go  to  1000 

else  if  (h  . eq .  2)  then 

x  1  0  =  x  1  3 
x  20  =  x  23 
y  0  =  y  3 
go  to  200 
end  i  f 

call  RSM  (xl3,x23,y3,xlf.x2f,yf) 

1000  print  *  ,  'the  end' 

end 

subroutine  TW0K(xl0,x20.y0,bl ,b2) 
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integer  j 

real  xl0,x20,y0,rf ,xln,x2n,xlp,x2p 

real  #y  1  1 ,*ylh,*yhh,*yhl ,bl , b2 , bl2 , bO , norm, m 

real  x 1 f , x2f , y f 

r  f  =0 . 5 

print  *,  ’Entering  2-k  factorial.’ 

print  »,  'The  radius  of  the  factorial  design  is’ ,rf 
print  *,  ’Would  you  like  to  change  this  radius  7 ’ 
print  *,  ’  1)  yes  2)  no’ 

read  * ,  j 
if  (j.eq.l)  then 
print  * ,  'Enter  the  new  radius 
read  *  . r  f 
end  l  f 

xln-xlO-rf 

x2n=x20-rf 

xlp=xlO+rf 

x2p=x20+rf 

call  S I M ( x 1 n , x2n , »y  1  1  ) 
print  »,  ’ » y 1 1  =’,»yll 

wr 1 1 e ( 7 , * )  , x In , x2n  ,  »y  1 1 

cal  1  S I M ( x 1 n ,x2p , «ylh) 
print  *,  ’»ylh=’,*ylh 
wr l t  e ( 7 , * )  ,xln,x2p,*ylh 

call  S I M ( x 1 p , x2n , *yh  1  ) 
print  *,  ’  »  y  h  1  = ’  . *  y h  1 

wr  l  t  e  ( 7  ,  »  ) ,xlp,x2n,*yhl 

call  S I M( x 1 p , x2p , »yhh ) 
print  *.  ’»yhh  = ’ , *yhh 
wr l te ( 7 , » ) ,xlp,x2p,*yhh 

m=  amax 1  ( »y 1 1  , *y 1 h  ,  *yh  1  ,  *y hh ) 
print  * ,  ’ m- ’ , m 

if  fyO.gt.m)  then 
print  *,  ’ yo  is  larger’ 

call  RSM  (xl0,x20,y0,xlf.x2f,yf) 
b  1  =0 
b  2  =  0 

go  to  70 
end  l  f 


b  0  =  !*yll+*ylh,-*yhl**yhh)  /  4 
bl  =  i  -  *yl  l-«ylh  +  «yhl  *-»yhh)  /  4 
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V' W'V'W 


'•»  '.■WL’*  V*'V"  V 


7”J»  V  'V  "V '  '.r  -V’  w :  v  „  w;  w\w  \  ■  •.•  *  v  ■  J"  *  -  «. ' 


* 


* 


<» 


70 


b2  = (-*yll+*ylh-*yhl 


b  12=  ( 
print 
print 
print 


*y  1 1  - *y 1 h - *yh 1 
* ,  ' b0= ’ , bO 

* ,  ' bl = ' , b 1 

* ,  ’ b2= ' , b2 


♦  *  y hh )  /  4 

♦  *  y hh ) / 4 


print  *,  ’ b 1 2  =  ’  , b 1 2 

wr Ite(7,*),’b0,bl,b2,bl2  are 
wr i te ( 7 , * )  ,  bO  ,  b  1  , b2 , b  1  2 


print  *,  ’note  interaction  of  bl2 


if  (  (bl  .eq.O)  .and.  (  b  2  .  e  q  .  0 )  )  then 
if  ( bO  . eq .  yO)  then 

print  *  ,  ’the  surface  is  flat  in  this  area' 
print  *  ,  ' y  = ’  ,  bO 

else  if  (yO  .gt.  bO)  then 
print  *  ,  ’max  point  in  this  area  is  '  . y  0 
else  if  (yO  .It.  bO )  then 
print  »,  yO ,  ’is  a  minimum  point  in  this  area’ 
end  i  f 

print  *  ,  'try  a  new  starting  point’ 
go  to  70 
end  if 


norm= (  ( b  1  *  *  2  )  +  ( b2  *  *  2 )  ) *  *  0  . 5 

b 1 :b 1 /norm 

b2=b2/norm 


continue 

return 

end 


subroutine  S I M ( x 1  , x  2 , y ) 
real  xl ,x2,y,w( 101 ,u,v 
integer  rep.i 
r  e  p  =  1 


c  pr int  *  ,  ’ En  ter  the  number  o  f  repe  titions  of  the 

s l mu  1  a  1 1 o  n  ’ 

c  print  *,  'wanted  at  this  point  to  reduce 

experimental  error, 
c  r e  ad  *  , rep 

y  =  0 

do  20  i=l, rep 
u=  x  1 
v  =  x  2 

c - this  is  the  place  to  insert  the  simulation. 

w (  i )  =100*  (  v  -  u  *  *  2  )  *  »  2  1*  (  1  - u )  **2 
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*  */*/  *.  “Va  **■  ***  •  *  •*-  a  •*.  ■*.  - 

I  *■%  Jn  ■VAvV  A  A  k  R aV  a  »  i.  mV  -  *-»*■*-«  A  -  .* 


l  u  Km  K  m.  K  M 


,  v\  , 

l  >  A  k  'll  fcV  a\.j>  A  kV  M*k 


the  next  line  changes  the  problem  to  a 
remove  the  c  in  column  1  of  line  18  to 
insure  the  c  13  in  column  1  of  line  18 
w  (  l  )  =  -  w  (  1  ) 


minimi zat ion 
minimize 
to  max l ml ze 


y  =  y  +  w  (  l  ) 
conti n  ue 
y  =  y  /  r  ep 

print  * , x 1 , x2 , y 

return 

end 

subroutine  LINE(xl0,x20,y0,bl  ,b2,zl  ,  z2  ,  y z  ,  l  > 
integer  i,n 

real  xl(0:  10)  ,  x2  (0  :  1  0 )  .  y (0 :  10)  ,bl,b2 
real  xl0,x20,y0,tl  ,t2.t3.c.d.e  ,  f  ,  z  1  ,  z  2 , y  z 
real  si  . s2  ,  s3 
n  =  1 


print  »,  'Entering  line  search 
1=0 

x 1  ( l )  =x  10 
x2 ( l )  =  x20 
y  (  i  )  =  y  0 


l  =  l  +  1 

xl(i)=xl(i-l)+(f2*»i)»bl) 
x2 ( l )  =x2 ( l -  1 )  *  (  (2**i)  *  b 2 ) 
s  1  =  x  1  (  l  ) 
s  2  =  x  2  (  l  ) 

call  SIM  I3l.32.s3! 
y  (  1  )  =  s  3 


print  *  ,  x  1  (  l 
print  * ,  x2 i l 
print  » ,  y(i) 
wr  1 1  e  l  7  ,  »  )  x  1 


(  l  )  .  x  2  (  i  )  .  y  (  l  ) 


if  <y(l)  .  It. v 0 )  then 
print  »  .  y 1  It  v  0 ’ 


i f ' n . eq . 1 ! then 
b l  =  -b  1 
b2=-b2 
1 3  =  y  (  1  ) 
n  =  2 

print  «.  ’ b 1  , b  2  ,  1 3  ,  n  are' 

print  *  ,  b 1  , b  2  ,  1 3  .  n 

go  to  5 

else  i f ( n . eq . 2) then 
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1 1  =  y  (  1 ) 

1 2  =  y  0 
c  =  -2 . 0 

print  *  ,  ’ t 1  , 1 2  ,  c  are' 

print  *  ,  t 1  , 1 2 , c 

go  to  15 
end  if 
end  if 

l  f  (  y (  l  )  .  g  t  .  v  (  l  -  1  )  1  then 

if  (l  . eq .  10)  then 

print  *  ,  ’the  surface  has  increased  for  10  steps' 

print  *  ,  'in  this  direction.  It  appears  to  go  to' 

print  *  ,  'infinity  for  a  opti ma 1  point.' 

print  »,  'start  with  a  new  point. ’ 

go  to  80 

else 

go  to  10 

end  l  f 
else 
i  =  i  +  l 

xl  (i) =xl  (  l  —  1 )  -  (  b  1  *  2  *  * (  i  -  2 )  ) 
x2  (  l  )  =x2  (  l -  1 )  -  ( b2  *  2  »  *  ( l -  2 )  ) 
call  SIM  (xl(l),x2(i),y(i)) 
print  *.  x 1 ( l ) 
print  *,  x  2 ( l ) 
print  *  ,  y <  l  ) 

write  1 7 , ♦)  x 1  ( l )  . x  2  (  1 1  .  y  <  l  > 
end  if 

if  (  y  (  l  )  . g  e .  y (  i -  2 )  >  then 

1 1  =  y  (  i  -  2  ) 

1 2  =  y  (  i  ) 

1 3  =  y  (  l  -  1  ) 
else 

t  1  =  y  (  i  -  3  ) 

1 2  =  y  (  i  -  2  ) 

1 3  =  y  (  l  ) 
end  if 

print  »  , t 1  . 1 2  ,  1 3 

c  =  2  »  *  f  l  -  2  ) 

15  print",  'hello' 

d  =  1 1  -  l  2  »  1 2  )  + 1 3 
l f  ( d . e  q  0 )  then 
d = . 000000 1 
end  l  f 

e=(-3.0»c»bl»tl)+(4.0»c»bl*t2)-(c*bl*t3) 
f  =  ( -3.0*c*b2»tl  )  +  <4  . 0 • c • b  2  »  t  2 )  -  (  c « b  2  »  1 3 ) 
print  *.c,d.e,f 
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V.  •C.'  .  *  .  -  .  '  •  'T-  . 


i  >  ■  *  *  *  »V»  W*  »s  > 


if  (  t  1  .  eq .  y ( 1 >  )  then 

z  1  -  x  1  (  i  )  -  0 . 5  »  e  /  d 
z  2  =  x 2  (  l  )  -  0 .  5*  f  /d 
else  if  (tl  .eq.  y(i-2) 
zl=xl  (i-2)  -  0 . 5  »  e d 

z 2  =  x 2 (  l  -2  )  -0 . 5*  f  /d 

else  if  (tl  .eq.  y  (  l  -  3 )  ) 
z  1  =  x  1  (  i  -  3 )  -  0.5*e/d 

z2-  x2(i-3)  -  0 . 5»f /d 

end  l  f 

call  SIM  lzl.z2.yz' 


then 


then 


0  print  *  , z 1  , z 2  ,  y z 

0  continue 

return 
end 

SUBROUTINE  PARTAN  (  x  1  0  .  x 2 0 . x  1  2  . x 2 2  . b 1  . f 2 > 

real  xl0.x20,xl2,x22.bl.b2.norm 

print  ».  'Entering  PARTAN.' 

b  1  =  x  1  2  -  x  1  0 

b2=  x22  -  x  2  0 

norm  =  (  ( b  1  »  * 2 )  ♦  (  b 2  *  *  2 )  )  *  * 0 . 5 

b  1  =  bl  'norm 
b2=  b  2 / n  o  r  m 

print  ».  'PARTAN  bl  ls’.bl 
print  *  ,  ' PARTAN  b2  ii  , b  2 

return 
end 

SUBROUTINE  FAIX 

x 1 0 . x  2  0 , x 1  1  , x  2 1  ,xl2,x22,bl  .b2,xl3,x23.y3) 

real  xl0,x20,xll  ,x21  , x 1 2 . x  2  2 . b 1  .b2.xi3.x2  3.y3 
real  mo , r 1  ,  r2  ,  r3  ,  r ,c , assu 

print  *  ,  'Entering  FAIX.' 

- - mo  is  the  slope  between  point  0  and  point  2  in 

pace 

mo  =  (x22-x20)  / ( x 1 2  - x 1 0 ) 

- is  the  distance  between  point  0  and  point  ! 

rl-sqrU  (xll-xl0)*»2  +  (x21-x20)»*2) 

- r2  13  the  distance  bet  we  e n  point  1  and  point  2 

r2  =  sqrt (  ' x 1 2-x 1  1 1  * *2  +  <x22-x21)»*2) 

- r  is  the  ratio  of  these  t  wo  distances 

r  =  r  1  /  r  2 

- r3  13  the  distance  between  point  0  and  point  2 
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V  '%P  1  *w  -Jf  -w  -  V  ■  i  » v 


TT 


v  ■  u  p  w  ■  L-"  ^  v  \r»  V""  v»  tr*  w"v  '.-v  >jv  •jwv  rw*rc,*'r_  v  *r 


r 3  =  sqr  t ( 

(  x22- 

x  2  0 )  **2+  ( x  1  2  ~  x  1  0 )  *  *  2 ) 

print 

# 

'  mo  = 

*  ,  mo 

print 

# 

’  r  1  = 

'  ,  r  1 

print 

# 

'  r2  = 

’  .  r  2 

print 

#  # 

'  r3  = 

'  ,  r  3 

print 

* 

’  r  =  ' 

,  r 

wr l te ( 7  , 

*)  ' r 1 ,r2 , r3 , r , mo  are  ’ 

wr  l  t  e  (  7  , 

»  )  r  1 

,  r  2  ,  r  3  ,  r  ,  mo 

c - c  1S  a  parameter  describing  the  eccentricity 

c  =  (  r  »  mo  +  1  )  /  (  r  »  mo  -  mo  »  »  2  ) 

c - assu  times  r3  is  the  assumed  acceleration  length 


assu=  (c*  (  (c  -  1  )  «  *  2  )  »  mo  »  *  2  ) / (  (  1  *  ( c  *  »  2 )  *mo*»2)  **2) 

print  *  ,  '  c  =  ’  ,  c 

write (7 , * )  ’ c= ' ,c 

wr lte(7,»)  ’assu=' .assu 

print  »,'assu=',assu 

xl3=xl2>assu» (  x  1  2  -  x  1  0  ) 

x23=x22+assu» (  x  2  2  -  x  2  0  ) 

call  SIM  (x 13 , x23 , y3) 

ret  urn 

end 


SUBROUTINE  RSM  (xl0.x20.ymm.xlf.x2f.yf) 
real  x  1  0 , x20 , x 1 p , x 1 n , x2p , x2n , r f 

real  *  y 1 1  ,ylm,  »  y  1  h  ,  y  m  1  ,y  mm ,  *  y  mh  ,  *  y  h  1  ,  y  h  m  ,  »  y  h  h 
real  bO , b 1  , b2 , b  1  1  , b22 . b  1  2 
real  x 1 f  , x2  f  , y  f 

r  f  =  0 . 5 

print  *,  ’Entering  3-k  RSM.' 

xlp=xlO+rf 
x ln  =  x  10-rf 
x2p=x20+rf 
x2n=x20-rf 

call  SIM  ( x In , x2n  ,  *yl  1  ) 
write(7.«)  x In , x2n  ,  *y 1  1 

call  SIM  ! x In , x20 , y 1m! 
wr  1  t  e  (  7  ,  #  )  x  1  n  ,  x  2  0  .  y  1  m 

call  SIM  (xln.x2p,*ylh) 
wr ite(7,»)  xln,x2p,»ylh 

call  SIM  (xl0.x2n.vml) 
wr l t  e ( 7 . * )  xl0.x2n.yml 

wr 1 1  e ! 7 , » )  x 1 0 , x  2  0 . y mm 
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call  SIM  (xl0,x2p,*y  mh ) 
wr  1 te i 7 , * )  xlO,x2p,»y  mh 

call  SIM  ( x  1  p  ,  x2n , *y h 1 ) 
wr  1  te (7 , » )  xlp,x2n,«yhl 

call  SIM  ( x  1  p , x20 , yhm) 
wr  ite (7 , * )  xlp,x2  0,yhm 

call  SIM  (x  1  p  ,  x2p , *yhh ) 
write (7,*)  xlp,x2p,*yhh 

30  bO=(*yll+ylm+*ylh+yml +  ymm+ *ymh+»yhl+yhm+*yhh)  /  9 . 0 

b 1 = (*yhl+yhm+«yhh-*yll-ylm-»ylh) /6.0 
b2=  ( *  y 1 h  + *  y mh  +  *yhh-*yll-yml-»yhl)  / 6 . 0 
bl  1= (»yll+ylm+»ylh  +  *yhl+yhm+*yhh- 
2  *  (  y  m  1  +  y  mm  +  »  y  mh  )  )  /  6 . 0 

b22=  (*yll+yml+*yhl+»ylh+*y  mh  +  »  y  h  h  - 
2  *  ( y 1 m  +  y mm  +  y hm)  )  /  6 . 0 

b  1  2  =  (*yll+*yhh-»ylh-*yhl)  /  4  .0 


print 

’  bO  =  ’ 

.  bO 

print 

*  , 

’  b  1  -  ’ 

,  b  1 

print 

*  , 

’  b2  =  ’ 

,b2 

print 

’  bl  1  = 

’  .  b  1  1 

print 

#  , 

’  b22  = 

'  .  b22 

print 

’  b  1  2  = 

'  ,b!2 

wr  l  te (7  ,  * )  ’  bO  ,  b  1  , b2  ,  b 1 1  ,  b22 , b 1 2  =’ 

wn  te (7 , *)  b0,bi,b2,bll,b22,bI2 

x  2  f  =  x  2  0  ♦  (  ( -b 1 »b  1  2 )  +  ( 2»b  1  1 *b2 )  )  / (  ( b  1  2 *b 1 2 )  - 

(  4  *  b  1  1  *  b  2  2  )  ) 

x  1  f  =  x  1  0  +  (  (  -  b  2  »  b  1  2  )  +  (  2  *  b  1  1  »  b  1  )  )  /  (  (  b  1  2  *  b  1  2  >  - 

(  4  *b  1  1  «b22 )  ) 

cal  1  SIM  (xlf.x2f.yf) 
print  *,  xlf,x2f,yf 

wr  l te  (7  ,  * )  ,  ’ the  final  point  l s '  . x 1 f  , x  2  f  , y  f 

if  ( 4  »  b  1  1  *  b  2  2  .It.  b  1  2  *  b 1 2 )  then 
print  *,  ’  ( x 1 f  ,  x  2  f )  is  a  saddle  point.' 

wr  l  te  (  7 , * )  ( x 1 f  .  x2  f )  is  a  saddle  point.' 

else  l f (  (bll.lt.O)  ,and(b22.1t.0))then 
print  *,  ( x 1 f  .  x  2  f )  is  a  max l mum  point, 

wr  i  t  e  (  7  ,  » )  ’  ( x 1 f  ,  x  2  f )  is  a  maximum  point.’ 

else  i  f  (  fbll.gt.O)  .and.  ( b  2  2 . g  t . 0  >  > then 
print  *  ,  '  ( x 1 f  , x  2  f )  is  a  minimum  point.’ 

wr  i te ( 7 , * )  '  ( x 1 f  , x2  f )  is  a  mini  mum  point, 

end  i  f 
return 
end 
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Appendix  B:  Program  Results 

Problem-  y=10-5!xl-3'»*2  -15(x2-|-'7>«»2 


point  0  =  10.00000  10.00000  -  4  5  7  0 . 0  0 

9.500000  9.500000  -4285.000 

9.500000  10.50000  -4795.000 

10.50000  9.500000  -4355.000 

10.50000  10.50000  -4865.000 

bO  ,  b 1  ,  b2 , b 1 2  are 

-4575.000  -35.00000  -255.0000  0.0C0C000E* 

slope  in  xl  direction  is  -0.1359800 
slope  in  x  2  direction  is  -0.9907116 

point  1  =  9.864020  9.009289  -4070.034 

the  new  slope  for  xl  is  0.9907116 
the  new  slope  for  x2  is  -0.1359800 

11.84544  8.737329  -4096.162 

7.882597  9.281249  -4085.385 

point  2  is  9.606606  9.044616  -4069.682 

partan  slopes  are  -0.38C7506  -0.9246778 

r 1  , r  2 . r  3 , r , mo  are 

0.9999996  0.2598276  1.033208  3.846704 

2 . 428566 
c=  3.000038 
assu=  2 . 4  198353E-02 

faix  point  3  is  9.597086  9.021497  -4057.933 

8.845104  7.195260  -3183.407 

7  322102  3.496549  -1736.066 

4.276097  -3.900873  -142.2109 

-1815912  -18.69572  -2157,813 

1.230093  -11.29830  -282  7932 

par  tan /line  point  3  is  2.999984  -6.999997 

10.00000 

point  3  is  2.999984  -6  999997  10.00000 

2.499984  -7.499997  4.999965 

2.499984  -6.999997  8.749922 

2.499984  -6.499997  4.999879 

2.999984  -7.499997  6.250043 

2.999984  -6.999997  10.^0000 

2.999984  -6.499997  6.249957 

3.499984  -7.499997  5.000122 

3.499984  -6.999997  8  750079 

3  499984  -6.499997  5  000036 

bO  . b 1  . b2  , b 1  1  , b22  ,  b  1  2  = 

6.666667  7.8837  074E-0  5  -4.2915344E-05  -  1  .  7  5  C  O  0  O 

-3 . 750000 

0 . 0000000E  +  00 

the  final  point  is  2.999995  -7.000003  10  000 

( x 1 f  . x  2  f  )  is  a  max i mum  po int . 
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Problem: 


100  t  x2  -  x  1 


x  1  )  - 


point  0=  5.000000 

4.500000  4.500000 

4 . 500000  5 . 500000 

5 . 500000  4 . 500000 

5 . 500000  5  500000 

bO . b 1  , b2 . b  1  2  are 
-4354^  50  -20254 . 00 

slope  in  xl  direction  is  -0. 
slope  in  x2  direction  is  9. 

3 . 009922  5 .  193968 

-0.9702346  5.596906 

1.019844  5 . 397937 

point  1  =  2.035566  5 

the  new  slope  for  xl  is  -9.9 
the  new  slope  for  x2  is  -0.9 
1 . 836598  3 . 306307 

1  .  43866  1  -0 . 6738498 

1  .  637629  1.3  16228 

point  2  is  1.853121 
partan  slopes  are  -0.8995146 
r  1  .  r  2  ,  r  3  .  r  .  mo  are 

2  .  9792  1  3  1  . 833909 

0  .  4856958 

c  =  3.234417 

assu=  0.3 167595 
faix  point  3  is  0.8563176 
point  3  is  0.8563176 
point  0=  0.8563176 

0.3563176  2.487431 

0.3563176  3.487431 

1  . 3563  18  2 . 48743  1 

1  .  3563  18  3 . 48743  1 

bO . b 1  . b2 . b  1  2  are 
-500.2607  343.3817 

slope  in  xl  direction  is 
slope  in  x2  direction  is  - 


.000000  -40 0 1 6  ' 0 

-24818.50 
-2  1768 . 50 
-66326  50 
-6 1276 . 50 

2025.000  -.0.- 

. 9Q5u39 . 

. 9484265E -02 
-1494.510 
-2171 . 297 

-  1899 . 092 

5.296385  -133.9797 

9484265E-02 

995039  1 

-  1  .  1459  1  8 
-752.9236 

-  186 . 8933 

3.471574  -0.8685583 

6  -0.4368905 


1 . 833909 


4984  1  9 


6  2  4  5  1  6 


48743  1 
487431 
487431 


2.987431  -508 

2.987431  -508.1403 

98743  1  -508  .1403 

-557 . 5954 
-1129. 689 
-42 . 09572 
-27  1.  6623 


343 . 38  17 
direction  is 
direction  is 


-200 
0 . 8636594 
-0 . 5040758 


4  15  1 


67  17- 


2  .  583637 
-0.8710013 
point  1  =  0 . 

the  new  slope 
the  new  slope 
-0 . 597  16  18 
-2.613465 
-  1  .  6053  1  3 
point  2  is  -0 
partan  slopes 
rl  ,r2,r3,r.mo 


1  .  979279 
3 . 995582 
4109897 
for  xl  is  -0 
for  x2  is  -0 
1 . 520028 
-  1  .9346  10 
-0 . 20729 1 1 
. 6546982 


-2207 . 654 
-  1051.  278 
3 . 247347 
5040758 
8636594 

-137. 9068 
-7695 .244 
-782.0326 
1 .421448 


-948.0226 


7068 


6943644 


196  2  36 
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••••  v 

DwUWllA  u‘iA  *_*.*.! V*Vy  *. . 


0.5156290 
. 036377 
c  =  -  1.525327 

assu  =  -0 . 8533999 

f  a  i  x  point  3  is 


2.114142 


0 . 6348025 


point  3  is 
point  0-  0. 

0.6148025 
0.6148025 
0 . 6548025 
0 . 6548025 
bO  , b 1  . b2 . b  1  2 
-554.5975 
slope  in  xl 
slope  in  x2 
2  .  206698 
-0 . 9370932 
point  1  =  0 


0.6348025 
6348025  2 

2 . 737858 
2.777858 
2 . 737858 
2.777858 

are 

11.97165 
direction  is  0 
direction  is  -0 
1 .521272 
3994443 
5173576 


the  new  slope 
the  new  slope 
-0.7192279 
-3  .  192399 
-1.955813 
point  2  is  -0 
par  tan  slopes 
-  1.8  16807 
0 . 9486414 


for  x 1  is  -  0 
for  x  2  is  -  0 
1  278354 
-  1  .  865437 
-0.2935417 
4340829 

are  -0.6913621 
0 .  19580  17 
3 . 085835 


2  .  176114 


2.757858 
2 . 757858 
757858 

-557 . 0496 
-576 . 0886 
-533 . 3094 
-  5  5  1  .942  1 


-9.4 17908 
. 7859479 
. 6 182927 

-1122.531 
-974 . 8847 

2.850250  -66' 

6  18  2  9  2  7 

7  8  5  9  4  7  9 

-60.877 7 9 

-  1  4554  .  33 

-  1705 .145 
1 .640818 

-0.7225083 
-972 . 0287 

-  477 .8248 


-  554 
-554 .6810 
554 .6810 


-  2 1 3  0  0  0  4 


partan/line  point  3  is  -0.  1003690 
393 . 0498 

point  3  is  -0.1003690  1.989566 

point  0=  -0.1003690  1.989566 


1  989566 

-393  0490 
-393.0498 


-0.1  1C3690  1  .  979566 

-0.1  103690  1  999566 

-9 . 0369038E-02  1.979566 

-9 . 0369038E-02  1.999566 

bO , b 1 . b2 . b 1 2  are 
-393.0207  -0.7726593 

slope  in  xl  direction  is  -0 
slope  in  x2  direction  is  -0 
-0.4834923  2.6605010E- 

-  1  249739  -3 . 8993  1  8 

-0  .  8666  157  -  1  .  936356 


-388 

-396 

-389 

-397 


2933 
2028 
8  306 
7562 


-3 . 958778 
. 1915617 
. 98 14806 
02  -6.492270 

-2987 , 493 
-725.6849 


•  e,G 


point  1 
the  new 
the  new  slope 
1  .  537095 
-2.398827 
point  2  is  0 
par  tan  slopes 
0  5005429 
-0 .  13040  14 


=  -0.4258662 
slope  for  xl 


0 . 32 18570 
is  0.98 14806 
for  x  2  is  -0.19156  17 
-6  .  1266333E-02  -587.8308 

07049803  -2513.000 

1850708  0  2026166 


4  on n. q 7  a, 


-  3  4  99  8  0  2 


are 


0  157736  1 

772346 
1 77579 


-0 

409 

468 


98748  1  3 
4  5  76 
086  1 


0 .  i 36  1  1  27 


partan/ line  point  3  is  0 
1 . 6037 18 

point  3  is  0.1956938 
-0.3043062  -0.3638873 

-0 . 3043062  0 . 1361127 

-0.3043062  0.6361127 

0.1956938  0 . 3638873 

0 . 1956938  0 . 1361127 

0 . 1956938  0.6361127 

0.6956938  -0.3638873 

0 . 6956938  0 . 1361127 

0.6956938  0.6361127 

bO . b 1  , b2  .  b  1  1  , b22  ,  b  1  2  = 

-2  1. 89624  -5  15  1964 

-25.00000 
19. 56938 


1956938 

0.1361127  -1.603718 

-22  .  53949 
-1.890530 

-  3 1  .24157 
-16. 82205 

-  1  . 6037  1 8 
-36.38538 
-7  1  .982  18 
-12. 19446 
-2.406738 

6.885005  -5  4387 


the  final  point  is  0.6847318  -2.5018558E-02  -24 

! x 1 f  . x  2  f  )  is  a  max l mum  point . 
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received  a  commission  in  the  USAF  through  the  ROTC  program. 
He  began  active  duty  in  November,  1974  at  Laughlin  AFB, 

Texas  in  the  Undergraduate  Pilot  Training  program.  Having 
completed  his  pilot  training  and  earned  his  wings  in 
‘’etcher,  1975,  he  served  as  a  KC-135  co-pi  lot  and  aircraft 
commander  at  Ellsworth  AFB,  S.  Dakota  until  May,  1981.  He 
attended  SOS  in  residence  in  summer  1981  and  was  assigned  to 
t he  4950th  Test  Wing  at  Wright-*Patterson  AFB,  Ohio.  His 
duties  consisted  of  research  pilot,  instructor  pilot  for  EC  - 
135  aircraft,  and  Chief,  Wing  Training  branch  until  entering 
the  School  of  Engineering,  Air  Force  Institute  of  Technology 
in  September,  1986 
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